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Dynamics of regular clusters of many non-touching particles falling under gravity in a viscous fluid 
at low Reynolds number are analysed within the point-particle model. Evolution of two families 
of particle configurations is determined: 2 or 4 regular horizontal polygons (called ‘rings’) centered 
above or below each other. Two rings fall together and periodically oscillate. Four rings usually 
separate from each other with chaotic scattering. For hundreds of thousands of initial configurations, 
a map of the cluster lifetime is evaluated, where the long-lasting clusters are centered around periodic 
solutions for the relative motions, and surrounded by regions of the chaotic scattering, in a similar 
way as it was observed by Janosi et al. (1997) for three particles only. These findings suggest to 
consider the existence of periodic orbits as a possible physical mechanism of the existence of unstable 
clusters of particles falling under gravity in a viscous fluid. 


I. INTRODUCTION 

In various contexts of the low-Reynolds-number fluid 
mechanics, periodic motions have turned out to be essen¬ 
tial for the dynamics of particle systems. Jeffery orbits of 
elongated particles or pairs of spherical particles in shear 
flow provide a classical example [TJ. There has been also 
a lot of interest in oscillatory motions of several parti¬ 
cles under external force fields. A modern example is a 
system of several particles kept by external forces inside 
a toroidal optical trap. Periodic oscillations, hydrody¬ 
namic pairing of particles, limit cycles, and transition 
to chaos have been found in experiments and theoretical 
models H-®. 

For many decades, periodic settling of a small number 
of particles under gravity has been extensively studied 
experimentally Emu and theoretically [13H2S], based on 
the Stokes equations for the fluid flow. A new perspective 
for the role of such benchmark solutions was opened by 
Janosi et al. EZj who showed that three particles settling 
under gravity in a viscous fluid in a vertical plane per¬ 
form chaotic scattering, and related this behaviour to ex¬ 
istence of an unstable periodic solution (which, however, 
was not directly found). For a wide range of arbitrary 
random initial configurations, the particle stay together 
if they are sufficiently close to this periodic solution; the 
closer they are, the longer they stay together m- Later, 
it was shown that three spheres exhibit a similar chaotic 
scattering [28j, and for spherical particles, the periodic 
solutions related to chaotic scattering have been explic¬ 
itly found [25] . 

It is important to check if a similar behaviour can be 
also observed in case of clusters made of a very large 
number of particles: is chaotic scattering of many par¬ 
ticles at non-regular configurations also observed, and is 
it coupled to periodic oscillations of regular arrays? Is 
the destabilization pattern of a suspension drop mm 
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related to the existence of periodic solutions for relative 
motions of the particles in certain regular configurations, 
as suggested in Ref. [52]? To address these open chal¬ 
lenging problems, the first step is to find and analyze 
examples of periodic oscillations of a large number of 
particles settling under gravity in a viscous fluid. 

In Ref. [16], periodic oscillations of a cube with the 
walls parallel and perpendicular to gravity were reported, 
and it was supposed that other regular polyhedrons may 
display similar behaviour. However, it has turned out 
that the periodic orbits of the cube m generalize to 
another family of periodic oscillations of an arbitrary 
large (even) number of particles: the particles arranged 
in two mirror regular horizontal polygons, one above the 
other [32] ■ For a similar geometry of arrays of rods (and 
other non-spherical particles), centered at vertices of a 
regular horizontal polygon, periodic solutions have been 
also found, both experimentally and theoretically ; 33] , 

In this paper, we start from generalizing the results of 
[32] . where periodic solutions were explicitly found for a 
moderate number of particles centered at vertices of two 
regular horizontal polygons (‘2 rings’). In this work, we 
explicitly demonstrate that 2 rings made of a very large 
number of particles indeed fall oscillating periodically. 
Next, we modify the initial particle positions to obtain 
‘less regular’ configuration, or in other words - we desyn¬ 
chronize the motion of the particles, and investigate if 
periodic solutions still can be found, and if they are re¬ 
lated to the existence of clusters with long lifetimes and 
chaotic scattering of particles in configurations which are 
close to the periodic orbits. 

The modification of the initial conditions is based on 
the idea of Ref. [32], where the 2-rings configurations 
were desynchronized by moving every second particle 
from each ring to the position it would have after one 
fourth of the period. Such a perturbation resulted in the 
initial configuration of four regular horizontal polygons 
(‘4 rings’), for which long-lasting quasiperiodic relative 
motions of particles were reported [32] . The question is 
if there exist some periodic solutions for an initial config¬ 
uration close by. Therefore, in this work we investigate 
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in details the dynamics of 4 rings for a wide range of the 
initial conditions and different numbers of particles. We 
evaluate a map of the cluster lifetimes, search for periodic 
solutions and chaotic scattering. 


II. SYSTEM AND ITS THEORETICAL 
DESCRIPTION 


We investigate dynamics of regular groups of many 


point-particles settling under identical gravitational 
forces G in a fluid of viscosity rj at the low-Reynolds- 
number. The fluid velocity v and pressure p satisfy the 
Stokes equations, see e.g. [33]. 
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where M is the number of particles, fi denotes the posi¬ 
tion of the particle with label i, and the 2 -axis is chosen 
along gravity, with G = — Gz where G > 0 and z is the 
unit vector along the 2 -axis. The equations are written 
in the laboratory frame of reference. 

The fluid velocity v(r) and pressure p(r) are given by 
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with the Green tensors, 
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where R = \R\. 

The equations of motion of the particles 

are given by 

M 
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where f, tJ = r, — fj, and u 0 = G/(6nija) is the Stokes 
velocity of the isolated particle. While solving Eqs. 0 , it 
is convenient to follow [323 anc i choose the inertial frame 
of reference moving with the velocity Uo- The benefit 
is that in this frame the dynamics is independent of the 
additional length scale a. The transformation back to the 
laboratory frame is straightforward, but not needed in 
this work, where we analyze the dynamics in the centre- 
of-mass frame of reference. 

From now on, we use dimensionless variables, based 
on an initial size of the group d as the length unit, and 
G/(8TTr]d) as the velocity unit. Therefore, 8mid 2 /G is the 


time unit. From now on, r, = fi/d denotes the dimen¬ 
sionless position of a particle i. 

Initially, the particles are placed at vertexes of K hor¬ 
izontal regular polygons (called ’rings’) which are sepa¬ 
rated from each other vertically and centered one above 
(or below) the other. The diameter d of the top ring is 
taken as the length unit. The number N of particles in 
every ring is the same. Therefore, the total number of 
particles in the system is equal to M = KN. Because of 
the system symmetry, we use the cylindrical coordinate 
system, in which r,, is represented as r, = (p it fa, 2 ,). 

Searching for periodic solutions, we take into account 
that numerical non-symmetrical perturbations can de¬ 
stroy periodic unstable solutions after times smaller than 
the period, as observed in Ref. [32). Therefore, we 
symmetrize the dynamics: we force the azimuthal com¬ 
ponents of the particle velocities to vanish, with <j>i = 
const(t) for i = 1 ,...,M, and the polygons to remain 
regular for all times. We use the parametrization i = 
K(n — 1) + k, with n = 1,..., N and k = 1,..., K. Then, 

27r(n-l) ^ 

<PK(n-X)+k ~ -jy-1 

and the radial and vertical coordinates of the particles 
from the same polygon k are the same, 


PK(n-V)+k — Pk, (9) 
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for n = and k = 1 The system of 3 M 

equations 0 is reduced to the system of 2 K equations 
for pi and zi , with l = 1 which in the frame of 

reference moving with the velocity Uo can be explicitly 
written as, 
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The numerical integration of Eqs. (11)-(13) is based 


on two methods: the fourth-order adaptive Runge-Kutta 
for non-stiff problems and backward differentiation for¬ 
mula (BDF) for stiff problems. We solve a system of 
ordinary differential equations using lsoda from the FOR¬ 
TRAN library odepack. Initially we start with non-stiff 
solver. Depending on number of time steps and accuracy 
the solver switches between method appropriate for stiff 
or non-stiff problem. This procedure is a standard and 
commonly used method for solving differential equations. 
This solver is sufficiently stable (35.:. The numerical cal¬ 
culations were performed with double precision and the 
error per each time step is not greater than 10~ 12 . 
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FIG. 1: (Color online) The system of M particles which form 2 rings, (a) M = 256. The initial configuration as in Eq. (151, 
(b) M = 20000. In the centre-of-mass frame, the trajectory of each particle has the same shape (not drawn to scale). Different 
shapes correspond to the indicated values of C. 


We are interested in relative motion of the particles; 
therefore, we will trace their positions in the centre-of- 
mass frame of reference, located at the symmetry axis, 
r GM = (0,0,2cm)- In this frame, pi and <fa are the 
same as in the laboratory frame, and particle vertical 
coordinates are Zi = Zi — zcm- 

We perform simulations for two families of the systems: 
with two and four rings. The specific initial configura¬ 
tions and their evolution will be discussed in the next 
sections. 


whole simulation time t = 5000 (what corresponds to 
300-2700 periods), see our movie [T] in 15B] . 

In the centre-of-mass frame, two particles with the 
same angular coordinates move along the same trajec¬ 
tory. Its shape is the same for all the pairs, and it de¬ 
pends on C as shown in figure [l](6) for M = 20 000. We 
found similar families of shapes for M = 256, analogical 
to those reported by [32] for M = 16, and resembling 
periodic solutions for rigid rods, discussed by [331 and 
[37]. 


III. DYNAMICS OF 2 RINGS 


IV. DYNAMICS OF 4 RINGS 


We first consider a system of M = 2TV particles which 
are grouped in two horizontal rings, each containing N 
equally spaced particles. Initially, the rings are identical 
and placed one exactly above the other. The diameter of 
each ring is equal to 1 and the initial distance C between 
the rings is a parameter in our simulations, as shown in 
figure [lja) . The initial positions r, = (pi,4>i,Zi) of the 
particles i = 1,..., M are 


/1 27r(n — 1) C\ 

\2 ’ N ’ 2 ) ’ 

/1 27r(n — 1) C\ 

\2’ N ,_ ~2 ) ’ 
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with n = 1,..., N. We conducted around 100 simulations 
for M = 256 and 20000 particles and different initial val¬ 
ues of C £ [0.05, 2.5]. For these values of the parameters, 
we observe periodic motion of the particles during the 


A. Initial configurations 


In the second system, M = 47V particles are grouped 
in four horizontal rings, each consisting of N particles. 
Initially, the rings labeled 1 and 3 are placed at 2 = 
C/2 and — C/2, respectively, and the rings labeled 2 and 
4: at z = 0, as shown in figure [2ja). The radii of the 
rings 1,3,2 and 4 are equal to 1/2, 1/2, R 2 , and i? 4 , 
respectively, with i? 4 > R 2 . The angular coordinates of 
the particles from ring 1 and ring 3 are the same. The 
other rings are rotated by ir/N around the symmetry 
axis. Summarizing, the initial cylindrical coordinates of 
the particles are, 
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FIG. 2: (Color online) Initially, M particles form 4 rings, as specified in Eq. (|19|) . Here, M = 256. 
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FIG. 3: (Color online) Evolution of M =256 particles which form 4 rings with C=1.5 and R 2 , R 4 as indicated. Trajectories 
of four particles (each from a different ring) in the centre-of-mass reference frame until time t = 500. The initial positions are 
marked by dots, (a) ‘3+1’ type of decay. One of the rings separates from the cluster, (b) ‘2+2’ type of decay. The group 
breaks up into two pairs of oscillating rings - one of them falls faster than the other, (c) ‘2+1+1’ type of decay. The cluster 
breaks up into one pair of oscillating rings and two single rings, (d) ‘4’: lack of decay. All rings oscillate. 


T An—X 


r±n 


f 1 27T (n — 1 ) C\ 

V2’ N ) ’ 





(18) 

(19) 


where n=l,...,N and R 2 and R 4 are the parameters of 
the simulation. 

We performed around 400 000 simulations for clusters 
made of different numbers of particles M = 64, 256,1024 
at the initial configurations specified by Eqs. 
with C £ [0.05,2.5], and a wide range of R 2 and R 4 . In 
Secs. |IVB||rVDl we will analyze how the dynamics de¬ 
pends on i ? 2 and R 4 for M = 256 and C = 1.5. Then, 
in Sec. |IVF[ we will argue that these results are generic 
also for the other values of M and C. 


B. Basic features of the dynamics 


The dynamics of four rings of particles is more complex 
than the dynamics of two rings. Although the general 
pattern of oscillations combined with settling is kept, in 
case of four rings the motion in general is not periodic and 
we observe destabilisation and decay of the system, as in 
movie [2j see [36]. Majority of initial conditions lead to 
the system decay during first 1000 units of the simulation 
time, what corresponds to about 100 oscillations. 

If the group breaks up, usually one ring is left behind 
the other three or one ring falls faster than the rest of 
the group; we call it ‘3+1’ type of decay. If the system 
separates into two pairs of rings, we denote it ’ 2 + 2 ’ type 
of decay. We observe also ’2+1+1’ type of decay when 
two rings oscillate together and the others are separated. 
Examples of each decay type are presented in figure (3|a)- 
(c) and movies S' -c in |36j . 
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For certain initial configurations, the particles perform 
quasiperiodic motion over a very long time. It has been 
checked that for several thousands of values of R 2 and 
R 4 the cluster lifetime exceeds 100000. A typical shape 
of the corresponding trajectory Zi(pi) is shown in figure 
[3^d) and movie ii, see [36]. These findings may indicate 
that for certain values of R 2 and JS 4 , periodic solutions 
exist, and they will be searched for later in this paper. 

We would like to remind that all simulations are per¬ 
formed for symmetrized configurations and the particles 
by definition do not change the vertical planes they be¬ 
long to, as described in Sec. [TT} The symmetrisation is 
done on purpose to find periodic solutions and ensure 
that the system will not break up because of non sym¬ 
metrical numerical perturbations before a period is com¬ 
pleted. 


C. Decay of the cluster 

To describe the system dynamics, it is crucial to define 
the event of the cluster decay and its time, denoted as 
r. The intrinsic property of a cluster of particles settling 
under gravity is the existence of periodic relative motions 
which may be used as an indicator whether the particles 
move together and interact with each other. For this 
reason we introduce a criterion of a cluster decay based 
on presence or absence of oscillations between pairs of 
particles. In particular, oscillations of particles i and j 
which belong to the same group imply that the difference 
of their z-coordinates, A Zij = Zi — z-j. oscillates around 
zero. We recognise that the particles i, j interact and 
stay together if Ahas repeating roots. 

Therefore, we use the following definition of the cluster 
decay. For each pair of particles i, j we calculate the dif¬ 
ference Aof their z-coordinates as a function of time. 
If for any l, m the time interval between two consecutive 
roots, tA and ts > tA, of A zi m exceeds 1000, we classify 
it as the cluster decay and denote min/ jm tA as the cluster 
lifetime r. 

This criterion of decay works well, in contrast to at¬ 
tempts based on measuring only the relative vertical dis¬ 
tances between rings, because it is common that rings 
approach each other again, even though they were sep¬ 
arated by a very long vertical distance meanwhile, as 
observed by [2Sj for systems of three particles. 

To define the type of the cluster decay we apply the 
following procedure: we count the number of roots of 
Ain the time interval [r, min(r + 1000, T)] for each 
pair of particles ij. If the number of roots is smaller 
than the length of the time interval divided by 200 (for 
the interval [r, r + 1000] it is equal to 5), we define that 
the particles i,j interact with each other. Otherwise, we 
denote them as separated. The number L of interacting 
pairs indicates the type of the decay. L = 3 corresponds 
to ‘3+1’ type of the decay, L = 2 to ‘2+2’ type of the 
decay, L = 1 to ‘2+1+1’ type of the decay and L = 0 to 
T+l+l+r. For L = 6 there is no decay (type ‘4’). 


Now, the main question is how the lifetime depends 
on initial conditions. The results are presented in fig¬ 
ure [4] as a map in space of the parameters R 2 and 1 ? 4 , 
with colours indicating the logarithm of the cluster life¬ 
time. Values of (R 2 , R 4 ) which lead to long-living clusters 
form a few compact regions, around which only isolated 
configurations with long lifetimes can be found. Deter¬ 
ministic regions are visible when R 2 is small enough, or 
J ?4 is large enough, and the cluster splits up immediately 
without oscillations. In the first case, settling velocity of 
a very small ring 2 is much larger than velocities of the 
other rings, while in the second case, settling velocity of 
a very large ring 4 is much smaller than velocities of the 
other rings. 

The essential property of the map shown in figure [4] is 
that for a wide range of R 2 and R ^ 1 the cluster lifetime 
is very sensitive to initial conditions. Around the regions 
of long-living clusters (with very sharp edges) we find a 
variability of the cluster lifetime by orders of magnitude. 
This property is also well visible in figure [5] at the cross 
sections of the map, plotted in the linear scale. This 
behaviour is similar to the chaotic scattering reported 
by [27] for three point particles sedimenting in a vertical 
plane. 

In figure [ 6 ] we show the type of the cluster decay and 
labels of the particles which interact with each other at 
the decay time r, depending on R 2 and R 4 . Sensitivity 
to initial conditions is clearly visible at most of the map 
areas, except a few deterministic regions corresponding 
to a very long or very short lifetimes. The last ones can 
be easily understood by comparing the ring diameters, 
and taking into account that the smaller the ring, the 
larger is its velocity. 


D. Clusters with long lifetimes and quasiperiodic 
solutions 

We have found that the long lifespan of the cluster 
and quasiperiodic trajectories of particles are observed 
for the system with initial configurations grouped in a few 
regions (see figure]!]). We will now investigate trajectories 
from each of the specific long-lifetime regions, and search 
for periodic solutions. 

To this goal, we apply the Dynamic Time Warping 
(DTW) method, described e. g. by (38] and [39] • For 
each particle i, we define its reference trajectory Zf(/ 9 j) 
for times t± <t < £ 3 , where tk, k = 1,2, 3, are the first, 
second and third time instants when Zi(tk) =0. Then, 
for f 3 < t < 5000, we calculate the Euclidean distance 
di(t) between the particle position (Zj(f), Pi(t)) and its 
reference trajectory. We define the average deviation A 
of the trajectories as the average of di(t) over all the 
particles i and all times t. 

In figure [4] we surround three main long-lifetime is¬ 
lands in the R 2 -R 4 space by rectangular boxes of a fixed 
resolution, and evaluate for each pixel the average devia¬ 
tion A of the trajectories. The results are shown in figure 
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FIG. 4: (Color) The cluster lifetime for around 420 000 initial configurations with C = 1.5 and different values of R 2 and R 4 
(drawn to scale). The darkest red colour corresponds to the cluster lifetime larger than the simulation time t = 5000. Resolution 
of the map is 0.002 x 0.002 and M = 256. 




FIG. 5: Cross-sections (now in the linear scale) through the map of the cluster lifetime from figure [4] (a) R 2 = 0.3, (b) 
Ri = 0.87. Sensitivity to initial conditions is clearly visible. Inset: Scaling down the resolution of R 2 and Ra from 0.002 to 
10~ 4 , a self-similar structure is found. 
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FIG. 6: (Color) Sensitivity of patterns of the cluster decay to the initial parameters R2 and R4, for C = 1.5 and M = 256. 
White colour means no decay, (a) The types of the decay: into ‘3+1’, ‘2+2’, ‘2+1+1’ or T+l+1+1’ rings (see Sec. |iVCl for 
the exact definitions), (b) ‘3+1’ type of the decay: the colour indicates the label of the ring which separates from the others 
at time r. (c) ‘2+1+1’ type of the decay: the colours indicate the interacting pair of the rings at time r. (d) ‘2+2’ type of the 
decay: the colours indicate which rings form interacting pairs at time r. 


[TJ For each of the three boxes, there exists a smallest de¬ 
viation A m j n , corresponding to a very thin trajectory. 
For example, in figure [7](a), A m i n = 3 • 10 -4 is three hun¬ 
dred times smaller than deviation of the quasiperiodic 
trajectories shown in figure [3j(d). 


E. Three periodic solutions 

The minima visible in figures [7|a), (6),(c) correspond 
to three periodic trajectories in the centre-of-mass frame 
of the cluster, shown in figures |8j<z), (6),(c), respectively. 
We remind that owing to the symmetrization of the dy¬ 
namics, each particle from ring i has the same coordinates 
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FIG. 7: (Color) The average deviation A of particle trajectories for clusters with lifetimes r > 5000, C = 1.5 and M = 256. 
Resolution: (a): 0.001 x 0.001; (6) and (c): 0.002 x 0.002. The minima at the subplots (a),(b),(c) correspond to the solutions 
presented in Fig.8(a),(b),(c), respectively. 
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FIG. 8: (Color online) Periodic trajectories of four particles (each from a different ring) in the centre of mass frame of a cluster 
with M = 256 particles. Dots: the initial particle positions, with C = 1.5. Solutions in the subplots (a)-(c) correspond to the 
minima of A in Figs. [7](a)-(c), respectively. 
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Zi(t) and Pi(t). Therefore, it is sufficient to display Zi(pi) 
for particles i = 1, 2, 3,4: each from a different ring i, with 
particles 1, 3 moving in vertical plane </>=0 and particles 
2,4 moving in another vertical plane <f> = 4.i t/M. Here, 
M = 256 and (7=1.5. 

The new periodic solutions have the following proper¬ 
ties. 

Solution 1, shown in figure (8^a) and movie i* in i| 36] . 
is located in the middle of the biggest of the long-lifetime 
regions from figure [4j and at the minimum of the trajec¬ 
tories deviation A in Fig. [fla). All four particles move 
along the same trajectory Zjfpf), shifted in phase by T/4, 
with the period T = 11.7. 

Solution 2, shown in ffgure|8jb) and movie^ in [136], 
is located in the middle of a small long-lifetime region 
from figure [4j and at the minimum of the trajectories 
deviation A in Fig. [7|b). It corresponds to the initial 
radius R 2 of the ring 2 much smaller than the initial 
radius of the rings 1 and 3, i?i = i? 3 = 0.5, and with R 4 > 
0.5. The central particle 4 has its own tiny trajectory 
and is circulated by the other three, which move along 
a larger trajectory with the period T = 10.7, shifted in 
phase by T /3 with respect to each other. The period of 
particle 4 is equal to T/3. 

Solution 3, shown in figure j8jc) and movie [8J; in I 3B] , 
is located in the middle of the other small long-lifetime re¬ 
gion from figure[4j and at the minimum of the trajectories 
deviation A in Fig. 7|c). The corresponding initial con¬ 
figuration satisfies the relations, R 2 ,Ra > R\ =Rz = 0.5. 
Now particle 2 is the central one, and the motion is qual¬ 
itatively similar to solution 2, with the interchange of 
the particles 2 and 4, and much wider rings, what leads 
to a larger period, T = 15.6, because larger rings move 
slower. 


F. Two families of periodic solutions 

Periodic oscillations have been found for many values 
of C. In figure [9] we show how the shape of the periodic 
trajectory of type 1 depends on C. Periodic solutions 1, 

2 and 3 for C = 2 are shown in figure 10 in the Appendix. 
These examples illustrate that the basic properties of the 
periodic solutions, described in the previous section and 
shown in figure [8j are generic for a wide range of values 
of C. 

Initially, two particles are at the same horizontal plane, 
and the other two particles are exactly one above the 
other. Movies [8 Jd,c in J3d illustrate that for each peri¬ 
odic solution 2 or 3, such a configuration appears again 
at T/6, but with different values of C, R 2 and i? 4 , which 
correspond to solutions 3 or 2, respectively. This new 
configuration can be treated as another initial condition 
for the same periodic solution. Therefore, assuming that 
there exist periodic solutions for all sufficiently large val¬ 
ues of C, we predict the existence of twin solutions 2 and 

3 for each value of C which is large enough. 

Dynamics of 4 rings made of other numbers M = 



- £7=0.5, R 2 =0.3702, i? 4 =0.7015 

- £7 = 1.5, R 2 =0.308, R t =0.821 

- £7 = 2.0, R 2 =0.304, R t =0.837 


FIG. 9: (Color online) Periodic solutions of type 1 for different 
values of C and M = 256. 


64,1024 of the particles have been also evaluated, and 
similar families of periodic solutions have been found. 


V. DISCUSSION AND CONCLUDING 
REMARKS 


In this work, we searched for periodic relative motions 
of many particles which form regular clusters with aspect 
ratios of the order of one. We generalized in this way 
the previous solutions found in Ref. [32] for a smaller 
number of particles. In Ref. [35], the systems of 2 rings 
were shown to destabilize. Therefore, to determine accu¬ 
rately periodic motions, in this paper we solved the sym¬ 
metrized dynamics, given in Eqs. (11)-(13). Although we 


have not analyzed here the non-symmetric perturbations, 
but it is known from Ref. |32j that such perturbations 
destabilize the system after times which are large enough 
to observe the existence of periodic solutions and deter¬ 
mine the cluster lifetimes (e. g. for 64 particles, periodic 
motions destablize at times t not shorter than half of the 
period, T/ 2, while even the motion during Tj 4 would 
be sufficient to deduce the existence of unstable periodic 
solutions). 

A similar behaviour was observed experimentally in 
Ref. [14j for three close particles sedimenting approxi¬ 
mately in a vertical plane. Typically, periodic motions 
have been observed during times comparable to one sixth 
of the period, T/6. Owing to symmetries of the motion, 
this is sufficient to demonstrate the existence of unstable 
periodic orbits. Moreover, in these experiments, the same 
type of periodic orbit was reached again after the destabi¬ 
lization, and this pattern of approaching and leaving the 
periodic orbit was repeated several times. Summarizing, 
examples of unstable periodic motions can be easily ob¬ 
served in experiments, and sometimes unstable periodic 
motions form a generic feature of the dynamics. 

In this work, we found a surprisingly large number of 
periodic solutions. For a cluster made of 2 rings, they 
were observed for all the investigated shapes (parame¬ 
terized by moderate values of the initial aspect ratio C ), 
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Solution 1: 
R 2 =0.304 
i? 4 =0.837 


Solution 2: 
R 2 =0.193 
i? 4 =0.547 


Solution 3: 
fi 2 =0.735 
i? 4 =1.148 


Z i 


Zi 


Zi 


particles 1 and 3, (f> 1 =</> 3 =0 particles 2 and 4, 4> 2 =4>4 =7r/64 



0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 

Pi Pi 


FIG. 10: (Color online) Periodic trajectories of four particles (each from a different ring) in the centre of mass frame of a cluster 
with M = 256 particles. Dots: the initial particles positions. Here, C = 2.0 (compare with figure [8] where C = 1.5). 


and particle numbers N even as large as 20 000. 

It has been interesting to modify the initial 2-rings con¬ 
figuration, with particles arranged in 4 rings, and in this 
way desynchronized. It turned out that clusters made of 
4 rings usually break up, and their lifetime and type of 
decay are in general sensitive to initial conditions. In the 
2D space of initial parameters, the deterministic area sur¬ 
rounds the chaotic region. The last one contains islands 
of quasiperiodic oscillations which do not destabilize dur¬ 
ing extremely long simulation times. At the centers of 
these islands, three different periodic solutions exist, pa¬ 
rameterized by the initial cluster height and the number 
of particles. 

The results provide a new perspective for the physical 
mechanism of decay of particle clusters sedimenting in 
a viscous fluid, typically with a very wide range of life 
times. 
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Appendix A: Periodic trajectories for a different 
value of the parameter C 

In this Appendix, we illustrate that Fig. [S] with 
C = 1.5 provides a generic example of two families of 
periodic solutions: solution 1 corresponds to the first 
family, and the solutions 2 and 3 to the second family, as 


discussed in Sec. IV F Such a structure of the periodic 


solutions has been observed for many values of the 
parameter C, which determines the initial geometry of 
the upper and lower rings. In Fig. |10[ we show another 
example: the periodic orbits for C = 2 (compare with 
Fig. [8|). 
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